function PV_cgns(grid,ubase,pbase,udofs,pdofs,up,k,outfile,test_name)
% PV_cgns(grid,ubase,pbase,udofs,pdofs,up,k,outfile,test_name)
%
% Write a (u,p) solution up of the Navier-Stokes problem (with
% continuous elements) in a format readable by paraview (.vtu).
%
% Velocity and pressure are interpolated on the same grid of degree k;
% k is independent from the degree of the finite element spaces.
% Bubble functions are handled correctly. The paraview output is
% structured as a discontinuous space, which simplifies the
% interpolations.
%
% Note: paraview only understands 3D data, and the usual way to work
% with 2D data is by setting the third coordinate to 0. This function
% handles the 2D and 3D case, adding zeros as required for the 2D
% case.
%
% See also: interpolate.

 % Interpolate pressure and velocity to the same grid
 for i=1:grid.d
   ui = up(i:grid.d:grid.d*udofs.ndofs);
   ui_disc = ui(udofs.dofs);
   [Xi,Ti,Ui] = el_resample(grid,ubase,ui_disc,k,1);
   U(i,:,:) = permute(Ui,[3 1 2]);
 end
 ui = up(grid.d*udofs.ndofs+1:grid.d*udofs.ndofs+pdofs.ndofs);
 ui_disc = ui(pdofs.dofs);
 [Xi,Ti,P,cell_type] = el_resample(grid,pbase,ui_disc,k,1);

 
 % Write the output file
 fid = fopen(outfile,'w');
 
 line = '<?xml version="1.0"?>';
 fprintf(fid,'%s\n',line);
 line = '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="BigEndian">';
 fprintf(fid,'%s\n',line);
 line = ' <UnstructuredGrid>';
 fprintf(fid,'%s\n',line);

 % Write the grid
 npoints = size(Xi,2)*size(Xi,3); % total number of points
 n_v4e   = size(Xi,2);    % number of vertexes per element
 n_v4c   = size(Ti,1);    % number of vertexes per cell
 n_c4e   = size(Ti,2);    % number of cells per element
 n_cells = n_c4e*grid.ne; % total number of cells
 line = [ '  <Piece NumberOfPoints="' , num2str(npoints,'%1d') , ...
                 '" NumberOfCells="'  , num2str(n_cells,'%1d') , '">'];
 fprintf(fid,'%s\n',line);

 line = '   <Points>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Float32" NumberOfComponents="3" Format="ascii">';
 fprintf(fid,'%s\n',line);
 if(grid.d==2) % add z to the data
   Xi(3,:,:) = 0;
 endif
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%f ',Xi(:,i,j));
     fprintf(fid,'\n');
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '   </Points>';
 fprintf(fid,'%s\n',line);

 line = '   <Cells>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="connectivity" Format="ascii">';
 fprintf(fid,'%s\n',line);
 for i=1:grid.ne
   i_shift = (i-1)*n_v4e;
   for j=1:n_c4e
     fprintf(fid,' %i',i_shift+Ti(:,j)-1);
     fprintf(fid,'\n');
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="offsets" Format="ascii">';
 fprintf(fid,'%s\n',line);
 i_shift = 0;
 for i=1:grid.ne
   for j=1:n_c4e
     i_shift = i_shift + n_v4c;
     fprintf(fid,' %i',i_shift);
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="types" Format="ascii">';
 fprintf(fid,'%s\n',line);
 for i=1:grid.ne
   for j=1:n_c4e
     fprintf(fid,' %i',cell_type);
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '   </Cells>';
 fprintf(fid,'%s\n',line);

 % Write the data
 line = '   <PointData Scalars="scalars">';
 fprintf(fid,'%s\n',line);

 % velocity
 line = '    <DataArray type="Float32" Name="velocity" NumberOfComponents="3" Format="ascii">';
 fprintf(fid,'%s\n',line);
 if(grid.d==2) % add Uz to the data
   U(3,:,:) = 0;
 endif
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%f ',U(:,i,j));
     fprintf(fid,'\n');
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);

 % pressure
 line = '    <DataArray type="Float32" Name="pressure" Format="ascii">';
 fprintf(fid,'%s\n',line);
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%f\n',P(i,j));
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);

 % vorticity
 switch grid.d
  case 2
   ui = cgns_vorticity(ubase,grid,udofs,up);
   ui_disc = ui(udofs.dofs);
   [Xi,Ti,O] = el_resample(grid,ubase,ui_disc,k,1);
   line = '    <DataArray type="Float32" Name="vorticity" Format="ascii">';
   fprintf(fid,'%s\n',line);
   for j=1:size(Xi,3)
     for i=1:size(Xi,2)
       fprintf(fid,'%f\n',O(i,j));
     end
   end
   line = '    </DataArray>';
   fprintf(fid,'%s\n',line);
 endswitch

 line = '   </PointData>';
 fprintf(fid,'%s\n',line);

 % Close the grid
 line = '  </Piece>';
 fprintf(fid,'%s\n',line);

 % Close the file
 line = ' </UnstructuredGrid>';
 fprintf(fid,'%s\n',line);
 line = '</VTKFile>';
 fprintf(fid,'%s\n',line);

 fclose(fid);
return

